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Abstract — The Discrete Fourier Transform (DFT) 
underpins the solution to many inverse problems com- 
monly possessing missing or un-measured frequency 
information. This incomplete coverage of Fourier space 
always produces systematic artefacts called Ghosts. In 
this paper, a fast and exact method for de-convolving 
cyclic artefacts caused by missing slices of the DFT is 
presented. The slices discussed here originate from the 
exact partitioning of DFT space, under the projective 
Discrete Radon Transform, called the Discrete Fourier 
Slice Theorem. The method has a computational com- 
plexity of 0(n\og 2 n) (where n — N 2 ) and is constructed 
from a new Finite Ghost theory. This theory is also 
shown to unify several aspects of work done on Ghosts 
over the past three decades. The paper concludes with 
a significant application to fast, exact, non-iterative im- 
age reconstruction from sets of discrete slices obtained 
for a limited range of projection angles. 

Index Terms — Discrete Radon Transform, Mojette 
Transform, Discrete Tomography, Image Reconstruc- 
tion, Discrete Fourier Slice Theorem, Ghosts, Number 
Theoretic Transform, Limited Angle, Finite Ghost The- 
ory 
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I. Introduction 

The Discrete Fourier Transform (DFT) is an important 
tool for inverse problems, where the Fourier representation 
of an object is used as a mechanism to recover that object. 
For example, in Computed Tomography (CT), the internal 
structure of an object is recovered or reconstructed from its 
projected "views" or projections 1 j. The Fourier Transform 
(FT) of the projections can be placed into Fourier space 
and the inverse FT is used to reconstruct the object. 
This is especially advantageous given the efficient and low 
computationally complex algorithm of the DFT known 
as the Fast Fourier Transform (FFT) [2 . In such cases, 
the acquired data cannot fully cover Fourier space as the 
problem is ill-posed [31 H]. Thus, there are missing Fourier 
coefficients and it is common practice to interpolate the 
space from the known Fourier coefficients. The choice of 
interpolation method is a major factor in determining the 
quality of the reconstruction [5] [6] . 

A. Ghosts 

Incomplete Fourier coverage leads to the introduction of 
reconstruction artefacts known as "Ghosts" or "invisible 
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Figure 1. An example of discrete and finite Ghosts, (a) shows a 
Ghost, as constructed by Katz [5], which is invisible when projected 
at any of the four rational angles shown, (b) shows a finite (periodic) 
Ghost, which is invisible in three finite angles (m = 1, 2, 3) shown. 



distributions" [7 j. Ghosts are effectively formed from an 
under-determined set of projections, i.e. from the "missing" 
or unmeasured projections. Therefore, Ghosts are always 
present within CT reconstructions and are the main reason 
for the filtering and interpolation of projection data [3J. The 
initial work on Ghosts in continuous space was pioneered by 
Bracewell and Roberts J7] , Logan [3] , Katz [8] and Louis [9] . 
Recent work by Candes et al [10 has also addressed this 
issue using more modern methods. 

In contrast, projection sets in Discrete Tomography 
(DT), which have practically well defined reconstruction 
processes [11], are often deliberately or naturally under- 
determined for various applications. These applications 
range from image encoding [12], network transmission [13] 
to tomography [T3] , Discrete Ghosts were first proposed 
by Katz [8 J as a way to describe zero-valued discrete 
projections taken at rational angles 6 pq , i.e. pq — tan _1 (g/p) 
(or simply the vector [q,p]) where p, q e Z. A simple example 
of these discrete Ghosts is shown in Fig. [TJa). Katz [8] 
determined that an N x N image can be reconstructed 
exactly from a set of p, rational angle projections if and 
only if 



N ^ 1 + max | ^2 \P: 

J=0 



/x-l 



(1) 



This is now known as the Katz Criterion. It is a state- 
ment that the information contained in the projection 
set needs to be one-to-one with the image data. Thus, 
knowing whether the projections of a sub-region in the 
reconstruction meets the Katz Criterion will allow one to 
determine if the sub-region is exactly reconstructable. 

Chandra et al. p~3J [15] showed that Ghosts also exist in 
the DFT as cyclic artefacts when DT projection data is 
missing. These cyclic Ghosts will be referred to as "Finite 
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Figure 2. An illustration of finite Ghosts formed within the Discrete 
Radon Transform (DRT). (a) shows missing projections (black rows 
in DRT space) and their effect on the reconstructed image, (b) shows 
the result of "De-Ghosting" in order to restore an image when a 
redundant image area is present. 




Figure 3. The Fourier eigenvalues of one dimensional circulant Ghost 
operators for N = 479. 



Ghosts" (see Fig. [TJb)). They showed it was possible to 
recover missing projection data exactly in DT given a 
sufficient redundant region in the reconstruction. However, 
the method could only solve for a very small number of 
missing projections while still being computationally viable. 
Their method was also based on empirical observations of 
the Finite Ghosts. 

In this paper, we present a fast method of recovering 
missing slices from DFT space that may be applied 
efficiently to very large sets of missing projection data. 
The method allows exact recovery of these missing slices 
by removing Ghosts formed within the DFT, a process 
that will be referred to as "De-Ghosting". The slices 
are equivalent to the projections of the Discrete Radon 
Transform (DRT). A schematic of this process is shown 
in Fig. [2] using the DRT. The method will utilise a new 
theory of Ghosts in the DFT constructed in Sec. [TTJ A 
consequence of the theory are the special one dimensional 
(ID) Ghost operators shown in Fig. [3] from which the 
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Figure 4. The slices within the geometry of the Discrete Fourier 
Transform (DFT) for a 5 x 5 image. Each colour represents a slice of 
a different slope with the DC coefficient centred (black). Note that 
each vector shown is computed modulo N. (a)-(e) shows the slices 
with slopes ^ m < 4 (mod 5) in DFT space, (f) shows the row sum 
(perpendicular) slice in DFT space. 



De-Ghosting is computed. The end result is a method 
to reconstruct an object exactly from a limited number 
of projections (with incomplete DFT coverage) without 
interpolation (see Sec. III). The incomplete DFT coverage 
in this case is restricted to the set of discrete slices of the 
DFT as described in the next section. It is hoped that 
a deep understanding of such methods will allow more 
efficient approaches to inverse problems that utilise the 
DFT. 



B. Discrete Fourier Slice Theorem 

The Discrete Radon Transform (DRT) provides an exact 
partitioning of DFT space in the form of finite or cyclic 
slices [16, TIP]. An example of the slice partitioning is given 
in Fig. |4j The DRT projections have the same geometry 
as the slices, i.e. they are taken as lines or congruences of 
the form 

y = rax + t (mod N), (2) 

for aniVxiV image where m is the slope, t is the translate 
and m,t,y,x e No- The lines in this geometry wrap around 
the image in both rows and columns. 

In this paper, the term "image" is used interchangeably 
with a two dimensional (2D) "object" and assume that N 
is prime, which defines the simplest form of the DRTQ In 
general, one requires a total of N + n /p projections for an 



^-See Hsung et < 
dyadic DRT. 



[T7] or Chandra et al. 18 for a discussion of the 
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Figure 5. The lines within the geometry of the Discrete Fourier 
Transform with slope 1 of a 5 x 5 image, (a) shows the wrapping 
behaviour of the lines in an image lattice. Each grey scale shows pixels 
of the same translate, (b) shows the equivalent circulant matrix for 
the lines. 



exact reconstruction for an N x N image, where N — p n 
(which includes powers of two). Thus, the projections are 
taken along the vectors [l,m] (with < m < N) and [0, 1], 
and the slices are placed along the vectors [— m, 1] and [1, 0] 
in DFT space (see examples shown as arrows in Fig. 

Chandra et al. [19j [20] showed that these projections 
(and lines) are equivalent to circulant matrices (or simply 
circulants) for all translates t (see Fig. [5|. 

Definition 1 (Circulant [21]). A circulant is an N x N 
matrix containing a unique row f(j) with j = 0, . . . , N — 1 
replicated on each row, but where each row is cyclically 
shifted (modiV) by a certain number of elements to the 
right. 

The property of circulants important here is that a 
circulant is diagonalised by the 2D DFT. Eq. ([3| shows 
an example of this property, where the values Xj are the 
eigenvalues of the circulant. 





(b) 20 Circulants (c) 40 Circulants 




(d) 60 Circulants (e) 80 Circulants (f) 100 Circulants 

Figure 6. Illustration of the Circulant Back-Projection of a 101 x 101 
image of Lena using different number of circulants. (a)-(f) show 
the effect of an increasing number of slices placed into Discrete 
Fourier Transform space on the image. The artefacts on the partial 
reconstructions are Finite Ghosts, which are the topic of Sec.[n] Using 
all (102) projections results in a perfect reconstruction of Lena. 
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Since each slice in the DFT is obtained by diagonalising its 
corresponding circulant using the DFT (see Eq. (|3j), the 
placing of each slice in DFT space adds its corresponding 
circulant to image space until the image is recovered. This 
is known as Circulant Back-Projection (CBP) [20] and an 
example of CBP is given in Fig. [6] When all projections are 
diagonalised, the resulting slices fully tile all of DFT space 
precisely once (see Fig. [4| allowing the exact reconstruction 
of the object. This is made possible by the prime-sized space 
of the DFT/image, since the gcd(ra, N) = 1 always. 

The DRT was independently developed by 
Grigoryan [22 , Bolker [23 , Gertner [24 and Fill 25 . The 
partitioning allows for the exact 0(n log 2 n) (with n — N 2 ) 
recovery of images from their projections and is known as 
the Discrete Fourier Slice Theorem (FST) [22j [21 \W\ . Do 
and Vetterli [27J TIP] utilised the Wavelet Transform on 
the DRT projections for image compression and de-noising. 
Chandra [20] extended the DRT to the Number Theoretic 
Transform (NTT) using the circulant theory of the DRT 



Figure 7. An example of Minkowski Addition or Morphological 
Dilation 0. (a) and (b) shows 2PSEs based on the vectors [q,p]. 
(c) shows the construction of a parallelogram from the Minkowski 
addition of (a) with (b). 



and showed that the discrete FST still applies. This 
Numb 3r-Theoretic Radon Transform (NRT) will prove 
critica in the De-Ghosting process later on in this work. 

X|ie Npy methods will rely on initial work done by 
Normand et al. [13]. They generalised the Katz Criterion 
to arbitrary convex regions by arguing that a region is 
reconslructable if and only if the size of the Ghost becomes 
too large to fit within the image region. To show this, they 
studied the spatial extents of discrete Ghosts, which were 
determined using Mathematical Morphology. 



C. Ghost Morphology 

Ghost structures are constructed using dilations of 2- 
Point Structuring Elements (2PSEs). A discrete 2PSE is 
the rational vector [q,p] as shown in part (a) and (b) of 
Fig. [7] A dilation by a 2PSE is the addition of the 2PSE 
to every point in the current structure [2H]- Equivalently, 
the dilation of A by B is given by the Minkowski Addition 
of A and B as 



A B = {a + b : a e A, b G B}. 



(4) 



An example of a dilation is shown in Fig. [7] For constructing 
Ghost structures, each missing projection has a 2PSE 
given by their rational angle or vector [q,p], so that B 
m Eq. @ is a 2PSE. The dilations are then applied using 
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Figure 8. An example of Ghost convolution kernels as an extension 
of Fig. [7] (a) and (b) shows kernels based on the vectors [l,mj] 
for each of the missing projections, (c) shows the construction of a 
parallelogram from the 2D convolution of (a) with (b). 



all 2PSEs of the missing projection angles to construct a 
Ghost structure. 

In the next section, the theory of Finite Ghosts is 
constructed that describes the artefacts caused by missing 
slices within the DFT. The artefacts formed within the 
redundant region of a reconstruction can then be used to 
recover the missing slices. 



II. Finite Ghosts 
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Figure 9. Ghosts in image space when missing projections are present 
in the Discrete Radon Transform. The Q x P image (Lena) is of 100 x 
100 pixels embedded in a image space TV = 113. Assuming one missing 
projection a at m a = 2, then the table shows the circulant nature of 
the Ghost artefacts that become embedded over the reconstructed 
image values. 



The approach of Normand et al. [13 is only adequate 
when one is concerned with the spatial extent of the Ghosts, 
such as in defining the reconstruction criterion. When one is 
also concerned about the values of the Ghosts, the dilation 
operation must be coupled with a negation of the signs 
within the Ghost structures. The negation ensures that 
the Ghost is also zero-valued in the direction of the 2PSE 
(see Fig. [8| . This allows Ghost construction to become a 
convenient 2D convolution. In the case of the DFT, the 
resulting operation is a 2D cyclic convolution and the 
discrete Ghosts are cyclic (or finite) also [T5] . 

Proposition 1 (Finite Ghosts). Each missing projection a 
at [1, m a ] in the DRT, which corresponds to a missing slice in 
the DFT, forms artefacts superimposed on the reconstructed 
image in the form of a circulant. The unique row of the 
circulant is —a. These artefacts are called Finite Ghosts. 

Proof: When projections are missing, the CBP is 
incomplete since there must be N + n/ v projections for 
an exact reconstruction. Each projection is equivalent to a 
circulant in the reconstruction process. Thus, the remaining 
missing contributions must be a superposition of <A£ number 
of circulants with shifts raj, where j = 0, . . . , fA£ — 1 and 9t 
represents the number of missing projections. The Ghosts 
are negative-valued since they are missing contributions in 
the reconstruction. ■ 
These Ghosts can be seen in Fig. [6] as slices are placed 
into DFT space. A schematic of the circulant form of the 
these Ghosts is shown in Fig. [9j The circulants effectively 
represent a 2D convolution of the missing projection based 
on the vectors [l,ra] of the DRT. In the next few sections, 
a method for defining this convolution is made and then a 
method for de-convolving or De-Ghosting these Ghosts is 
constructed. 



A. Ghost Convolution 

Let the 2PSE of the DRT projections be defined as [1, mj] 
for each of the j projections. 

Proposition 2 (Ghost 2PSEs). The 2PSEs of the Finite 
Ghosts are [l,raj] for each of the missing projections j — 
0,...,fA£-l. 

Proof: From Prop, [I] each missing projection is a 
circulant shifted by raj and hence each artefact is shifted 
by raj on subsequent rows. Therefore, this is equivalent to 
the Minkowski Addition of each row by [l,raj]. ■ 

Proposition 3 (Ghost Kernels). The Ghost convolution 
kernels are the 2PSEs [l,raj] ; which sum to zero along its 
vector for each of the missing projections j = 0, . . . , *AC — 1. 

Proof: The projections of Ghosts are zero- valued at 
certain finite angles [l,raj] by construction (see Prop. [2|, 
so the 2PSEs must also be a Ghost at its finite angle, since 
they apply to each point in the image/structure. These 
then become Ghost convolution kernels (see Eq. Q). An 
example of Ghost kernels can be seen in Fig. [8j ■ 
Then the Ghost kernels can be used in an initial 
method to construct Finite Ghosts by computing the 2D 
convolution of the Ghost kernels as follows: 

1) Compute the DFT of the initial structure and each 
of the 2PSEs. 

2) Multiply the coefficients of each 2PSE, together with 
the coefficients of the structure, in DFT space to 
compute the 2D convolution. 

3) Compute the inverse DFT to get the Ghosts in image 
space. 

The method has a computational complexity of 
0(ftCN 2 log 2 N) for computing the 2D DFTs and 0(^N 2 ) 
for computing the convolution, where <A£ denotes the 
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number of Ghosts. An example of the 2D convolution 
method is given in Fig. [8| The resulting Ghost will be 
invisible for those projections at the finite angles used to 
construct them. 

However, this 2D convolution approach is inefficient 
because it does not take into account the following: 

1) The Ghost grows a row at a time and, in most cases, 
does not take up the full 2D space. Thus, a series of 
full 2D convolutions can be expensive. 

2) The DRT space contains zero valued projections and 
when 0>i is close to AT, this space is mostly empty. 

Is there a way to make the process more efficient? Yes, the 
answer is to use the Projection Convolution Theorem of 
the DRT. 

B. Projection Convolution Theorem 

The Projection Convolution Theorem (PCT) states that 
a 2D convolution is equivalent to the ID convolution of each 
projection or slice in Fourier space. Thus, the 2D Ghost 
convolution can be computed as a series of ID convolutions 
on the known projections. For Ghost convolution, this is 
particularly useful when the number of Ghosts fA£ is close 
to N, so that the number of known projections is small. 
Therefore, DRT space is sparse and the 2D convolution 
becomes a relatively small number of ID operations. 

The PCT can be interpreted as a consequence of the 
discrete FST. Since the slices tile DFT space exactly 
and that slices are the DFT of the projections, cyclically 
convolving the slices of two objects is equivalent to cyclically 
convolving the objects themselves. Thus, in order to utilise 
the PCT, one needs to know the projections of the two 
objects. In this case, the projections of the Ghost kernels 
are particularly convenient, which makes a PCT approach 
very efficient. 

Proposition 4 (Kernel Projections). The projections of 
the 2D convolution kernel [l,mg\, with a positive term at 
the origin and a negative term at the coordinate (l,mg), 
will have the positive term at the zeroth translate and the 
negative term at translate tj as 



(mg 



(mod AT) 



(5) 



for each projection mj with j = 0, . . . , N — 1 and Q = 
0, . . . , N — 1. For the j = N projection, negative term is 
always at tj = 1. Note that negative values —x modulo N 
are equivalent to the value N — x (mod N) . 

Proof: The projection mg will be zero by definition 
with all other projections being non-zero. For projections 
with mj < mg, the negative term will appear in the 
translates t > (modiV) because the initial 2PSE [l,raj] 
for t = will not sample the point (l,mg) while projecting. 
The difference of slope (mg — mj) defines the translate 
where the point is eventually sampled. For projections with 
mj > mg, the negative term will appear in the translates 
t < (mod Af) because the initial 2PSE [l,m] will not 
sample the point (l,mg) while projecting until the 2PSE 
wraps around the image. Hence, the difference of slope 




(a) 

Figure 10. A schematic of the ID Ghost convolution operators (a). 
Each row represents the unique projection of a Ghost convolution 
kernel, (b) shows their eigenvalues in (centred) discrete Fourier space. 



(mg 



will be negative and defines the translate where 



the point is eventually sampled. ■ 

Corollary 1 (Kernel Operators). The projections of the 
kernels [l,mg] only have N distinct combinations of positive 
and negative values. 

Proof: This follows from Prop. [4] and the fact that the 
system is constructed within a finite geometry, i.e. because 
of the wrapping of the values since there are only N residue 
classes. ■ 
Thus, the unique kernel projections (following Cor. [T]) 
and their eigenvalues have the form as shown in part (a) of 
Fig. [lO] These eigenvalues can be precomputed, since they 
depend purely on the set of finite angles of the missing 
projections as given by Prop. [4] The form of the DFT 
eigenvalues are shown in part (b) of Fig. 10 and Fig. [3] for 
N = 479. 

Thus, the ID Ghost convolution approach may be used 
to generate a 2D Ghost as follows: 

1) Pre-compute the Ghost operator eigenvalues via the 
DFT. This can be used as a hash table to pick out 
the relevant operator based on the projection being 
convolved. 

Convolve a delta function with each of the 2PSEs 
by selecting the correct eigenvalues for each operator 
using the eigenvalues hash table (as given by previous 
step) via Eq. ([5| and multiplying these operators with 
the eigenvalues of the delta function. 
Inverse DFT to obtain the Ghost structure in image 
space. 

The computational complexity is 0(fA£/Vlog 2 iV) for pre- 
computing the eigenvalues, 0(/nN \og 2 N) for the ID DFTs 
and 0(p9^N) for computing the convolutions. This method 
is advantageous when /i<A^, since the computations are 
done only on the small number of known projections. A 



2) 



3) 



simple example of this construction is shown in Fig. 11 for 
two 2PSEs. 

However, both methods are susceptible to round-off 
errors and loss of precision for large fA£ when utilising the 
DFT for the convolutions. This problem, which manifested 
as numerical overflow, was also encountered by Chandra et 
al. Q3] with their method to remove Finite Ghosts. This 
numerical growth can be easily seen as a direct consequence 
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(h) 



(i) 



Figure 11. A simple example of the ID convolution approach to 
constructing 2D Ghosts, where grey and black denote zero and minus 
greyscale values respectively, (a) and (d) show the 2PSEs [1, 1] and 
[1,2] respectively, (c), (f) and (i) show the DRT projections of the 
2PSE (a), the 2PSE (d) and the convolution of 2PSEs (g) respectively. 



of the convolutions of the eigenvalues in Fig. [3] The solution 
is to use the Number Theoretic Transform and the Number- 
Theoretic Radon Transform of Chandra [20] . 

C. Number Theoretic Convolution 

The Number Theoretic Transform (NTT) allows one to 
compute convolutions just as the DFT because the unit 
circle is replaced with the digital "circle" 

M-l 



1 (modM), (6) 

M_1 1 is a multiple of M, N is a multiple of M — 1 
[29] . The successive powers {1, . . . , M — 1} of 



so that a 
and a, M G 

a generates a unique set of integers in some order modulo 
M. Such a number a is called the primitive root of M. 

The primitive root(s) a in these cases have to be found 
by trial and error and can be computed by dividing 
M — 1 by the prime factors pj of M — 1, such that 
a (M-i)/pj ^ ^ (modM), where the trial value of a is prime. 
Integer coefficients allow computations to be done without 
round-off error or numerical overflow, since the results 
are congruent modulo M [30 . A new efficient and fast 
algorithm for computing prime-length NTTs, that may be 
utilised for this work, is described in [A] 

Chandra [20^ showed that the discrete FST still holds 
within the 2D NTT when placing Number Theoretic slices, 
i.e. the NTTs of the projections, into 2D NTT space. 
This allows one to replace the DFT in all computations, 
including those within the Ghost convolutions, with the 
NTT. The resulting NRT was constructed specifically to 
remedy the loss of precision when forming Finite Ghosts. 
Chandra [20] also showed that the implementation of the 




Figure 12. The Number Theoretic Transform eigenvalues of one 
dimensional circulant Ghost operators for N = 479. Note that the 
fine structure present is more easily seen when the above image is 
viewed from a distance. 





Figure 13. Examples of Ghosts in the Number-Theoretic Radon 
Transform, (a) shows fA£ = 4 within a TV = 13 with Ghosts at m = 
1, 2, 3, 4 (mod 53). (b) shows 9i = 80 within a TV = 101 with Ghosts 
at m = 1, . . . , 80 (mod 607). Each Ghost occupies fA£ + 1 rows in the 
image. 



NTT is faster than the DFT because of its integer-only 
operations. 

Consequently, the Ghost convolution method is im- 
pervious to numerical overflow and faster than DFT 
approaches. The NTT integer eigenvalues of the Ghosts 

The Ghosts in 



kernel projections are shown in Fig. 12 
this space have the value (mod M) in the direction of the 
missing projections. An example of Ghosts within the NRT 
can be seen in Fig. [l3j Ghosts in the NRT have another 
important property in that any physical or perceivable 
structures in the Ghosts are difficult to discern. This makes 
the NRT Ghosts well suited for encoding and encryption. 

In summary, the advantages of the ID Ghost convolution 
approach are 

1) Simplicity: ID operators are simple and cyclic convo- 
lutions are much easier to compute. 

2) Efficiency: The majority of DRT or NRT space will 
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Figure 14. To partially de-convolve the Q x P sub-image (in grey), 
one aligns the first row A of the Ghost constructed (in red) to the 
image row to be De-Ghosted. The remaining rows would be in the 
redundant region. Once (I) is completed, the De-Ghosted sub-image 
row is removed from the Ghosts in (II) to make it into a redundant 
row. The process is continued until all image rows are De-Ghosted. 



A visual interpretation to constructing Ghosts, that 
unifies the Ghost Recovery algorithm of Chandra [H] and 
the 2D convolution approach of this paper, can be made 
using n-gons. See the thesis of Chandra [31J for initial work 
on this topic. 

More work has to be done with both methods when 
DFT space has noise or inconsistencies present as these 
also become convolved during the De-Ghost process. The 
convolution approach requires a very good estimation of 
noise prior to De-Ghosting, so that the estimates may be 
used to de-convolve their effects on the results. Further 
work needs to be done in generalising the convolution 
approach to arbitrary missing discrete Fourier coefficients 
using the theory outlined in this work. Recent work by 
Svalbe et al. [32] discussed the minimal extent of Finite 
Ghosts, which may prove useful in this endeavour. In the 
next section, the De-Ghost method is applied to the discrete 
inverse problem of determining a reconstruction from a set 
of rational angle aperiodic projections. 



be empty due to the Ghosts, so one always keeps the 
amount of work to a minimum. 
3) Speed: ID blocks are at most N elements apart in 
memory. 

Once the Ghost is constructed/convolved, it can be 
used to remove the artefacts formed by the missing slices 
(represented via their 2PSEs) from which it is constructed. 
This process of de-convolution is described in the next 
section. 

D. Ghost De- Convolution 

The Ghost convolved in the previous section may be 
applied as a mask to remove the Ghost introduced by 
missing projections in the DRT or missing slices in the 
DFT. The mask is applied by individually convolving each 
of the rows of the Ghost with the rows of the image as 
shown in Fig. [141 



Part (I) of Fig. 14 is effectively a sum of ID convolutions 
or a sub-circulant expression given as 

A-/(3,0) T + B-/(4,0) T + C-/(5,0) T + J D-/(6,0) T = 7(3, 0) T , 

( ? ) 

where A,B,C and D are circulants, I denotes the Ghosted 
image and I denotes the desired De-Ghosted sub-image. 
Row 7(3, 0) is back-substituted into 7(3, 0) order to make 
it a redundant row. Then Eq. ^ is repeated for (II) and 
so on. In this method, it is important to note that the 
shifts need to be —mg, rather than mg, when convolving 
Ghosts to undo the right shift of the circulants. This is 
done by setting N - mg as the projection angle for the 
missing projections mg in the Ghost convolutions. 

The coefficients of the Ghosted image and Ghost struc- 
ture need only be computed once and the remaining 
computations done in DFT or NTT space. Thus, the 
computational complexity of this de-convolution method 
is O(QfAtp). Note also that the de-convolution requires a 
total of fA£+ 1 rows, with fA£ of those rows being redundant, 
to function. 



III. Discrete Reconstruction 

Chandra et al. Q3] utilised their Ghost removal technique 
to reconstruct an image from a limited angle set of 
discrete ID projections. The term limited angle refers to 
the fact that the coverage of the half-plane need not be 
uniform and/or complete. Recent work done by Chandra et 
al. [18] has been done in improving the use of the discrete 
projection data within the DFT for fast reconstruction. 
However, the angle set utilised in [18 covered the half- 
plane. In this section, the De-Ghost method will be used 
to extend the work of Chandra et al. [18] to include any 
limited angle set of projections and reducing the total 
number of unique projections (for an exact reconstruction) 
to Q + 1 for a Q x P image. 

For a Q x P image in a N x N space, a total of 
Q + 1 projections need to be taken because one needs 
fAC redundant rows in the image for fA£ missing slices. This 
means that there must be \i number of projections so that 

— N — /j, = N — (Q + 1) , since there are N + 1 total number 
of slices. The rational angle projection geometry is known 



as the Mojette Transform (MT). Fig. 15 shows a simple 
example of a MT for a 4 x 4 image using three projections. 



A. Ghost Angle Sets 

Chandra et al. [H] showed that the MT projection set 
directly and efficiently maps to a prime-sized DFT space 
exactly as 

m = pq~ X (mod TV), (8) 

where q~ x is the multiplicative inverse of q that can be 
computed easily via the Extended Euclidean algorithm. 
For the dyadic or power of two case, extra s projections 
are needed (see Chandra [20]), which are mapped as 



2s = p 1 q (mod N), 



(9) 



when the gcd(g, N) > 1, Positive [q,p] values represent 
the first octant of the half-plane with the other octants 
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Figure 15. An example of a Mojette Transform for a discrete image 
of size 4x4 using the three projections [1,1], [1,-1] and [1, —2]. The 
bold lines within the right-hand grid shows a possible reconstruction 
path using a corner-based reconstruction method |33) . 



produced by [— q,p], [p, q] and [p, —q] [53] . The reconstruction 
is then obtained with a computational complexity of 
0(n log 2 n) with n = iV 2 , which is the same order as the 2D 
FFT. 

The same mapping can be used to cover a fraction of the 
half-plane by limiting the number of octants used within the 
projection set. One then obtains a set that is limited, but 
covers DFT space fully because there exists a multiplicity 
of possible rational vectors for each finite angle [l,m]. This 
multiplicity is shown as graphs in Fig. [16] for the case 
of a quadrant and a half-plane. One may observe that 
the multiplicity is relatively "flat" in a discrete sense, but 

i is reminiscent of curves 



with small variations. Graph 16 



obtained by Svalbe and Kingston [35] when observing 
the "unevenness" of the rational vectors from uniform 
coverage. The unevenness of the Farey vectors is related 
to the distribution of prime numbers and the Riemann 
Hypothesis pi [37]. 

The theory of Finite Ghosts allows one to reduce the 
number of projections required to Q + 1, rather than a 
total N + 1, for these limited angle sets by choosing the 
projections with a desired property, such as the ones with 
the smallest number of bins. An example of such projection 




Multiplicity Histogram of the Finite Angles for N=256 in One Quadrant 




150 200 250 

Finite Angle 



(a) 



Multiplicity Histogram of the Finite Angles for N=256 in the Half-Plane 



50 100 150 200 250 300 350 400 

Finite Angle 

(b) 

Figure 16. The multiplicity of the mapping between the finite angle 
and rational angle sets for N = 256, i.e. for the dyadic case. Because 
of the dyadic size, the range of values are ^ m < N and s' = N + s, 
where ^ s < N/2. 



sets are given in Fig. 17 



B. Results 

Let us consider a 100 x 100 image of Lena and 101 rational 
angle projections similar to the geometry of Fig.[T7^b). Also, 
for performance, let N — 257, so that N — 1 is a power of 
two (see [A]) . Therefore, from the 101 projections, there will 
be 156 Ghosts to remove and also 156 redundant rows in 
the image, so that the De-Ghost method can be applied. 

Fig. [l8^a) shows the DRT space resulting from 101 
rational projections acquired across the half-plane and 
mapped via Eq. Fig. 18{b) and (c) show the resulting 
Ghosts superimposed on the reconstruction due to the 
missing projections in (a). Fig.[l8^d) shows the NTT Ghost 
eigenvalues that are used to construct the de- convolution 




(a) 



(b) 



Figure 17. Examples of rational projection sets when utilising Ghosts 
in discrete reconstruction, (a) shows a set limited to a quadrant and 
(b) shows a set using the half-plane. The sets apply to a 11 x 11 image 
embedded into a 23 x 23 space. 



Ghost in Fig. 18 'e). Fig. [18m shows the final result of 
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used to solve the discrete inverse problem of reconstructing 



(e) (f) 

Figure 18. The results of the De-Ghost method of Sec. [H] applied to 
discrete reconstruction. A total of 101 rational angle projections of 
a 100 x 100 image of Lena was exactly reconstructed within a DFT 
space of size 257 x 257. (a) shows the resulting DRT space, (b) shows 
the resulting Ghosts when (a) is reconstructed, (c) shows a cropped 
version of (b) focussing on the embedded image, (d) and (e) show 
the Ghost eigenvalues and the Ghost that was used to recover the 
reconstructed image exactly (shown in (f)). 



the De-Ghosting, which is an exact reconstruction. The 
computation time was approximately 500 milliseconds on a 
1.6 GHz AMD™ Turion™ 64-bit Laptop. This fairs a few 
orders of magnitude better than the method of Chandra et 
al. [E] , whose computations took in the order of hours to 
complete. 

Further work needs to be done to handle inconsistencies 
within projections while utilising De-Ghosting. Also, a 
study of the small variations within the rational angle 
multiplicity may be important in understanding whether 
some projections are more important than others. 

Conclusion 

A new discrete theory of Ghosts was constructed that 
allowed one to describe the exact effects of missing periodic 
slices in the Discrete Fourier Transform (DFT). The theory 
was used to construct new Ghost convolution and de- 
convolution methods, which can be utilised in recovering 
missing slices in the DFT (see Figs [8] to [l4| and Props [I] to [5]). 
The methods required the use of the Number Theoretic 
Transform and the Number-Theoretic Radon Transform of 
Chandra [20 to avoid numerical overflow and other loss 
of precision problems. This De-Ghost method was then 



from limited rational angle projections (see Figs 17 and 18). 
This solution has a computational complexity of 0(n log 2 n) 
with n — N 2 and results in an exact reconstruction. A brief 
study of the limited angle set was also made (see Graphs 



in Fig. 16). 
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Appendix 

For a prime-length iV, the NTT can be computed by first 
selecting the modulus M as M = k • N + 1. This allows the 
power of a in Eq. Q to be a multiple of N. For example, 
the modulus for the prime-length N = 101 is 607 with k = 6. 
Then Rader's [38] can be used to compute its fast form. 

Rader [38] constructed a Fast DFT algorithm for prime- 
length DFTs by turning them into cyclic convolutions 
of composite lengths. This is achieved by permuting the 
input vector x so that the Transform matrix F becomes 
a circulant matrix W. The algorithm is computed by 
handling the DC term separately and then applying the 
following four steps: 

1) Compute the primitive root a for the length N — p. 
Determine the new order x(j) — >> x(f) of the elements 
of the input vector where j' is given by j' = a _J 
(modp) = aP~i (modp). 

Determine the unique row of W the circulant form 
of V. This row is given by Wj = a aJ (modp). 
If p — 1 is not highly composite, determine a power 
of two N f > 2 • p - 4. Create a new sequence x' of this 
length and pad it with N 1 - p + 1 zeros between the 
first and second elements of the vector. Create W' 
of length N r by repeating the p — 1 sequence of W . 

An example of the circulant form W , as a 4 x 4 sub-matrix, 
is given in Eq. ( 10 ) for p — 5 and a ■ 
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The zeros in the above matrices are present because the 
DC term is computed separately and the first element of 
the input vector is added to the result after the convolution. 
The rest of the transform is completed as the following: 

1) Determine the eigenvalues of (or diagonalise) x and 
W' by computing the NTT of each using the dyadic 
divide and conquer algorithms [2j EH] • 

2) Multiply the eigenvalues to determine the eigenvalues 
of the convolution. 
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3) Inverse NTT the convolution eigenvalues to determine 
the convolution. 

4) Normalise the convolution by p — 1 using the multi- 
plicative inverse so that (p— 1) • (p— = 1 (modM). 

5) Add the first element of x to each coefficient of the 
transform result. 

6) Compute the DC term by summing all terms of the 
input vector and placing this as the u — coefficient 
of the transformed result x(u). 

7) Extract the rest of the coefficients of x(u) for < u < 
p — 1 of the prime-length transform elements by using 
u = a? (mod p) from the normalised convolution. 

8) Normalise the result by p if desired using multiplic- 
ative inverse so that p • p~ x = 1 (mod;/). 

The implementation of this transform can be found in the 
Finite Transform Library (FTL) library [30] . 
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